Sexual conversion rate - IgG and IgM reactivity measured by flow cytometry
2026-01-29
# load required libraries
library("lme4")
library("ggplot2")
library("here")
library("tidyverse")
library("emmeans")
library("brms")
library("DHARMa")
library("glmmTMB")
library("performance")
library("scales")1 Exploration and visualisation
A pseudonymized/codified dataset was used as a precaution due to ethical considerations (patient-level clinical data). The full dataset is available upon request (see manuscript data availability statement).
# load the data
df <- read_csv(here("./data/Ig_detection_codified.csv"), show_col_types = FALSE) %>%
# rename(Experiment = `Experiment (= Sample collection Number)`) %>%
rename(HBB.genotype = `HBB Genotype`,
`Raw IgM+ RBC counts_IgM plate` = `Raw IgG+ RBC counts_IgM plate`) %>%
mutate(across(c(`Date Flow cytometry experiment`), as.Date)) %>%
select(
ID, HBB.genotype, Sample_type,
`Raw IgG+ RBC counts_IgG plate`, `Raw IgG+ iRBC counts`, `Raw IgG+ Asexual iRBC counts`, `Raw IgG+ Sexual iRBC counts`,
`Raw IgM+ RBC counts_IgM plate`, `Raw IgM+ iRBC counts`, `Raw IgM+ Asexual iRBC counts`, `Raw IgM+ Sexual iRBC counts`,
`Raw IgG- iRBC counts`, `Raw IgG- Asexual iRBC counts`, `Raw IgG- Sexual iRBC counts`,
`Raw IgM- iRBC counts`, `Raw IgM- Asexual iRBC counts`, `Raw IgM- Sexual iRBC counts`,
Age_group, Gender, Village, CollectionSite,
Pf_infection, Parasites_µL, Gametocytes_µL, iRBCs_IgG_PercentageCells, iRBCs_IgM_PercentageCells
) %>%
rename(
raw_igGpos_RBC = `Raw IgG+ RBC counts_IgG plate`,
raw_igGpos_iRBC = `Raw IgG+ iRBC counts`,
raw_igGpos_asexual = `Raw IgG+ Asexual iRBC counts`,
raw_igGpos_sexual = `Raw IgG+ Sexual iRBC counts`,
raw_igMpos_RBC = `Raw IgM+ RBC counts_IgM plate`,
raw_igMpos_iRBC = `Raw IgM+ iRBC counts`,
raw_igMpos_asexual = `Raw IgM+ Asexual iRBC counts`,
raw_igMpos_sexual = `Raw IgM+ Sexual iRBC counts`,
raw_igGneg_iRBC = `Raw IgG- iRBC counts`,
raw_igGneg_asexual = `Raw IgG- Asexual iRBC counts`,
raw_igGneg_sexual = `Raw IgG- Sexual iRBC counts`,
raw_igMneg_iRBC = `Raw IgM- iRBC counts`,
raw_igMneg_asexual = `Raw IgM- Asexual iRBC counts`,
raw_igMneg_sexual = `Raw IgM- Sexual iRBC counts`,
) %>%
mutate(HBB.genotype = case_when(
HBB.genotype == 0 ~ "HbAA",
HBB.genotype == 1 ~ "HbAC",
HBB.genotype == 2 ~ "HbAS",
HBB.genotype == 3 ~ "HbCC",
HBB.genotype == 4 ~ "HbSC",
)) %>%
mutate(Pf_infection = case_when(
Pf_infection == 0 ~ "NO",
Pf_infection == 1 ~ "YES",
)) %>%
mutate(across(c(ID, HBB.genotype, Sample_type, Gender, Village,
CollectionSite, Pf_infection, Age_group),
as.factor)) %>%
mutate(HBB.genotype = fct_relevel(HBB.genotype, "HbAA")) %>% # set HbAA as baseline reference
drop_na(starts_with("raw_"))
all(df$HBB.genotype == df$Sample_type)1.1 IgG
plot_func <-
function(data, x, y, title){
data %>%
ggplot(
aes(x = {{x}},
y = {{y}}))+
geom_boxplot(aes(fill = {{x}}), outlier.size = 2, outlier.colour = "grey",
width = 0.8, position = position_dodge(width = 0.8)) + # Reduce space between boxplots
geom_point(aes(fill = {{x}}), alpha = 0.8, shape = 21, colour= "#115e67a6",
position = position_jitter(width = 0.1, seed = 42)
) +
scale_fill_manual(values =
c("HbAA" ="paleturquoise","AA" ="paleturquoise",
"HbAC" = "palegreen",
"HbAS" = "thistle",
"HbCC" = "lightpink1",
"HbSC"= "peachpuff")) +
scale_y_continuous(
# limits = c(0, 1),
labels = label_percent()) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 14), # Increase font size for x-axis ticks
axis.text.y = element_text(size = 14), # Increase font size for y-axis ticks
axis.title.x = element_text(size = 16), # Increase font size for x-axis title
axis.title.y = element_text(size = 16), # Increase font size for y-axis title
strip.text = element_text(size = 12),
panel.grid.major = element_blank(), # Remove major gridlines
panel.grid.minor = element_blank(), # Remove minor gridlines
axis.line = element_line(color = "black") # Add black axis lines
) +
labs(title = title)
}
plot_func(data = df, x = HBB.genotype, y = raw_igGpos_iRBC / (raw_igGpos_iRBC + raw_igGneg_iRBC), "iRBC IgG+ / iRBCs")
plot_func(data = df, x = HBB.genotype, y = raw_igGpos_asexual / (raw_igGpos_asexual + raw_igGneg_asexual), "Asexual IgG+ / asexual iRBCs")
plot_func(data = df, x = HBB.genotype, y = raw_igGpos_sexual / (raw_igGpos_sexual + raw_igGneg_sexual), "Sexual iRBC IgG+ / sexual iRBCs")
1.2 IgM
plot_func(data = df, x = HBB.genotype, y = raw_igMpos_iRBC / (raw_igMpos_iRBC + raw_igMneg_iRBC), "iRBC IgM+ / iRBCs")
plot_func(data = df, x = HBB.genotype, y = raw_igMpos_asexual / (raw_igMpos_asexual + raw_igMneg_asexual), "Asexual IgM+ / asexual iRBCs")
plot_func(data = df, x = HBB.genotype, y = raw_igMpos_sexual / (raw_igMpos_sexual + raw_igMneg_sexual), "Sexual iRBC IgM+ / sexual iRBCs")
1.3 Pf infection
ggplot(df,
aes(x = Pf_infection, y = raw_igGpos_iRBC / (raw_igGpos_iRBC + raw_igGneg_iRBC))) +
geom_boxplot(width = 0.8, position = position_dodge(width = 0.8)) +
geom_point(alpha = 0.8, shape = 21, colour= "#115e67a6",
position = position_jitter(width = 0.1, seed = 42))
ggplot(df,
aes(x = Pf_infection, y = iRBCs_IgG_PercentageCells/100)) +
geom_boxplot(width = 0.8, position = position_dodge(width = 0.8)) +
geom_point(alpha = 0.8, shape = 21, colour= "#115e67a6",
position = position_jitter(width = 0.1, seed = 42))
ggplot(df,
aes(x = Pf_infection, y = raw_igMpos_iRBC / (raw_igMpos_iRBC + raw_igMneg_iRBC))) +
geom_boxplot(width = 0.8, position = position_dodge(width = 0.8)) +
geom_point(alpha = 0.8, shape = 21, colour= "#115e67a6",
position = position_jitter(width = 0.1, seed = 42)) +
scale_y_continuous(
limits = c(0, .15),
labels = label_percent())
ggplot(df,
aes(x = Pf_infection, y = iRBCs_IgM_PercentageCells/100)) +
geom_boxplot(width = 0.8, position = position_dodge(width = 0.8)) +
geom_point(alpha = 0.8, shape = 21, colour= "#115e67a6",
position = position_jitter(width = 0.1, seed = 42)) +
scale_y_continuous(
limits = c(0, .15),
labels = label_percent())
Comparison between percentage and raw values.
plot_func(data = df, x = Pf_infection, y = raw_igGpos_iRBC / (raw_igGpos_iRBC + raw_igGneg_iRBC), "iRBC IgG+ / iRBCs")## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
plot_func(data = df, x = Pf_infection, y = iRBCs_IgG_PercentageCells, "iRBCs_IgG_PercentageCells")## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
plot_func(data = df, x = Pf_infection, y = raw_igMpos_iRBC / (raw_igMpos_iRBC + raw_igMneg_iRBC), "iRBC IgM+ / iRBCs")## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
plot_func(data = df, x = Pf_infection, y = iRBCs_IgM_PercentageCells, "iRBCs_IgM_PercentageCells")## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
Figure for publication:
# Color palette
fill_colors <- c("NO"="#FFDAB5",
"YES" ="#7BAFD4")
outline_colors <- c(
"NO"="#FF8C00",
"YES" ="#375D81")
# Updated outline colors (darker tones)
outline_colors <- c(
"NO"="#FF8C00",
"YES"= "#375D81"
)
summary_stats <- df %>%
group_by(Pf_infection) %>%
summarise(
ymin = min(raw_igGpos_iRBC / (raw_igGpos_iRBC + raw_igGneg_iRBC), na.rm = TRUE),
ymax = max(raw_igGpos_iRBC / (raw_igGpos_iRBC + raw_igGneg_iRBC), na.rm = TRUE),
ymed = median(raw_igGpos_iRBC / (raw_igGpos_iRBC + raw_igGneg_iRBC), na.rm = TRUE)
)
p <- ggplot(df, aes(x = Pf_infection , y = raw_igGpos_iRBC / (raw_igGpos_iRBC + raw_igGneg_iRBC) )) +
geom_jitter(aes(color = Pf_infection, fill = Pf_infection),
width = 0.15, size = 3.5, shape = 21, alpha = 0.9, stroke = 0.8) +
geom_errorbar(data = summary_stats, aes(ymin = ymin, ymax = ymax, x = Pf_infection),
width = 0.2, color = "grey60", inherit.aes = FALSE) +
geom_point(data = summary_stats, aes(y = ymed, x = Pf_infection),
shape = 95, size = 8, color = "black", inherit.aes = FALSE) +
scale_fill_manual(values = fill_colors) +
scale_color_manual(values = outline_colors) +
# scale_y_continuous(limits = c(0, 30), breaks = seq(0, 30, 10), expand = expansion(mult = c(0.03, 0.03))) +
scale_y_continuous(
# limits = c(0, .15),
labels = label_percent()
) +
geom_hline(yintercept = 0, linetype = "dotted", color = "grey40") +
theme_minimal(base_size = 15) +
theme(
panel.grid = element_blank(),
axis.line = element_line(color = "black"),
axis.text.y = element_text(color = "black", size = 14), # Bigger Y-axis numbers
axis.text.x = element_text(color = "black", size = 14), # Bigger X-axis category names
axis.title = element_text(color = "black", size = 16), # Bigger axis titles
legend.position = "none"
) +
labs(x = "Pf infection status", y = "Total IgG+ Troph-iRBCs (% cells)")
p
ggsave(here("./figures/pf_infection_IgG.png"))## Saving 7 x 5 in image# Summary stats
summary_stats <- df %>%
group_by(Pf_infection) %>%
summarise(
ymin = min(raw_igMpos_iRBC / (raw_igMpos_iRBC + raw_igMneg_iRBC), na.rm = TRUE),
ymax = max(raw_igMpos_iRBC / (raw_igMpos_iRBC + raw_igMneg_iRBC), na.rm = TRUE),
ymed = median(raw_igMpos_iRBC / (raw_igMpos_iRBC + raw_igMneg_iRBC), na.rm = TRUE)
)
# Color palette
fill_colors <- c("NO"="#FFDAB5",
"YES" ="#7BAFD4")
outline_colors <- c(
"NO"="#FF8C00",
"YES" ="#375D81")
# Updated outline colors (darker tones)
outline_colors <- c(
"NO"="#FF8C00",
"YES"= "#375D81"
)
p <- ggplot(df, aes(x = Pf_infection , y = raw_igMpos_iRBC / (raw_igMpos_iRBC + raw_igMneg_iRBC) )) +
geom_jitter(aes(color = Pf_infection, fill = Pf_infection),
width = 0.15, size = 3.5, shape = 21, alpha = 0.9, stroke = 0.8) +
geom_errorbar(data = summary_stats, aes(ymin = ymin, ymax = ymax, x = Pf_infection),
width = 0.2, color = "grey60", inherit.aes = FALSE) +
geom_point(data = summary_stats, aes(y = ymed, x = Pf_infection),
shape = 95, size = 8, color = "black", inherit.aes = FALSE) +
scale_fill_manual(values = fill_colors) +
scale_color_manual(values = outline_colors) +
# scale_y_continuous(expand = expansion(mult = c(0.03, 0.03))) +
scale_y_continuous(
# limits = c(0, .15),
labels = label_percent()
) +
geom_hline(yintercept = 0, linetype = "dotted", color = "grey40") +
theme_minimal(base_size = 15) +
theme(
panel.grid = element_blank(),
axis.line = element_line(color = "black"),
axis.text.y = element_text(color = "black", size = 14), # Bigger Y-axis numbers
axis.text.x = element_text(color = "black", size = 14), # Bigger X-axis category names
axis.title = element_text(color = "black", size = 16), # Bigger axis titles
legend.position = "none"
) +
labs(x = "Pf infection status", y = "Total IgM+ Troph-iRBCs (% cells)")
p
ggsave(here("./figures/pf_infection_IgM.png"))## Saving 7 x 5 in image2 Final models
2.1 IgG
2.1.1 IgG+ iRBC
Proportion of IgG positive versus IgG negative infected red blood cells.
A binomial GLM is not applicable because of violated model assumptions.
model <- glm(
cbind(raw_igGpos_iRBC,
raw_igGneg_iRBC) ~ HBB.genotype,
family = binomial, data = df)
summary(model)##
## Call:
## glm(formula = cbind(raw_igGpos_iRBC, raw_igGneg_iRBC) ~ HBB.genotype,
## family = binomial, data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.545394 0.009922 -155.762 < 2e-16 ***
## HBB.genotypeHbAC -0.020901 0.014188 -1.473 0.141
## HBB.genotypeHbAS -0.219481 0.014104 -15.562 < 2e-16 ***
## HBB.genotypeHbCC 0.207452 0.026210 7.915 2.47e-15 ***
## HBB.genotypeHbSC -0.012781 0.021102 -0.606 0.545
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 8242.3 on 50 degrees of freedom
## Residual deviance: 7807.8 on 46 degrees of freedom
## AIC: 8237.9
##
## Number of Fisher Scoring iterations: 4# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)
# Test for overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 8.0713, p-value < 2.2e-16
## alternative hypothesis: two.sided# Test for zero-inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)
The beta-binomial GLM was used as the final model. The data was modelled using generalized linear model using the beta-binomial family and a single dispersion parameter to account for overdispersion.
model <- glmmTMB(
cbind(raw_igGpos_iRBC,
raw_igGneg_iRBC) ~ HBB.genotype,
family = betabinomial(link = "logit"),
data = df
)
summary(model)## Family: betabinomial ( logit )
## Formula: cbind(raw_igGpos_iRBC, raw_igGneg_iRBC) ~ HBB.genotype
## Data: df
##
## AIC BIC logLik -2*log(L) df.resid
## 746.6 758.2 -367.3 734.6 45
##
##
## Dispersion parameter for betabinomial family (): 27.1
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.52094 0.12918 -11.774 <2e-16 ***
## HBB.genotypeHbAC -0.10308 0.18413 -0.560 0.576
## HBB.genotypeHbAS -0.24702 0.17918 -1.379 0.168
## HBB.genotypeHbCC 0.24564 0.34265 0.717 0.473
## HBB.genotypeHbSC -0.03087 0.27385 -0.113 0.910
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)
# Test for overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.81732, p-value = 0.42
## alternative hypothesis: two.sided# Test for zero-inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)
The estimated odds ratios (using Wald Z-tests) for the different genotype contrasts are as follows:
# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(model, ~ HBB.genotype)
# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## HbAC / HbAA 0.902 0.166 Inf 0.546 1.49 1 -0.560 0.9807
## HbAS / HbAA 0.781 0.140 Inf 0.479 1.27 1 -1.379 0.6414
## HbAS / HbAC 0.866 0.158 Inf 0.527 1.42 1 -0.791 0.9333
## HbCC / HbAA 1.278 0.438 Inf 0.502 3.26 1 0.717 0.9527
## HbCC / HbAC 1.417 0.488 Inf 0.554 3.62 1 1.013 0.8494
## HbCC / HbAS 1.637 0.559 Inf 0.645 4.16 1 1.442 0.6001
## HbSC / HbAA 0.970 0.266 Inf 0.459 2.05 1 -0.113 1.0000
## HbSC / HbAC 1.075 0.296 Inf 0.507 2.28 1 0.262 0.9990
## HbSC / HbAS 1.241 0.338 Inf 0.590 2.61 1 0.793 0.9326
## HbSC / HbCC 0.758 0.303 Inf 0.255 2.25 1 -0.692 0.9582
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 5 estimates
## Intervals are back-transformed from the log odds ratio scale
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)
contrast(emm, method = list("HbAC / HbAA" = AC - AA,
"HbAS / HbAA" = AS - AA,
"HbCC / HbAA" = CC - AA,
"HbSC / HbAA" = SC - AA),
adjust = "bonferroni",
type="response",
infer = c(TRUE, TRUE)
)## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## HbAC / HbAA 0.902 0.166 Inf 0.570 1.43 1 -0.560 1.0000
## HbAS / HbAA 0.781 0.140 Inf 0.499 1.22 1 -1.379 0.6720
## HbCC / HbAA 1.278 0.438 Inf 0.543 3.01 1 0.717 1.0000
## HbSC / HbAA 0.970 0.266 Inf 0.489 1.92 1 -0.113 1.0000
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 4 estimates
## Intervals are back-transformed from the log odds ratio scale
## P value adjustment: bonferroni method for 4 tests
## Tests are performed on the log odds ratio scaleThe above table not only shows the estimated odds ratio and its
corresponding p-value, but also the confidence interval of this odds
ratio (= asymp.LCL and asymp.UCL).
These p-values are adjusted for multiple testing using the Bonferroni method.
The estimates are already transformed from the log-scale, so they can be directly interpreted as odds ratios (OR).
The estimated SCR per group is, adjusted for multiple testing using the Bonferroni method:
confint(emm, type="response", adjust="bonferroni")2.1.1.1 Bayesian estimates
brm_model <- brm(
bf(
raw_igGpos_iRBC | trials(raw_igGpos_iRBC + raw_igGneg_iRBC) ~
HBB.genotype
),
family = beta_binomial(link = "logit"),
data = df,
iter = 5000,
chains = 4,
# control = list(adapt_delta = 0.95),
seed = 42,
cores = 8
)
## Compiling Stan program...
## Start samplingsummary(brm_model)## Family: beta_binomial
## Links: mu = logit
## Formula: raw_igGpos_iRBC | trials(raw_igGpos_iRBC + raw_igGneg_iRBC) ~ HBB.genotype
## Data: df (Number of observations: 51)
## Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
## total post-warmup draws = 10000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.52 0.14 -1.80 -1.25 1.00 9384 7216
## HBB.genotypeHbAC -0.10 0.20 -0.49 0.29 1.00 9731 7487
## HBB.genotypeHbAS -0.24 0.19 -0.62 0.14 1.00 9421 8149
## HBB.genotypeHbCC 0.18 0.40 -0.67 0.89 1.00 10567 6076
## HBB.genotypeHbSC -0.06 0.31 -0.69 0.50 1.00 11056 6669
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi 24.64 4.97 15.80 35.40 1.00 9980 7170
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(brm_model, ~ HBB.genotype)
# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))## contrast odds.ratio lower.HPD upper.HPD
## HbAC / HbAA 0.903 0.586 1.30
## HbAS / HbAA 0.783 0.499 1.10
## HbAS / HbAC 0.868 0.571 1.24
## HbCC / HbAA 1.230 0.417 2.27
## HbCC / HbAC 1.368 0.469 2.55
## HbCC / HbAS 1.578 0.563 2.94
## HbSC / HbAA 0.955 0.466 1.60
## HbSC / HbAC 1.057 0.509 1.78
## HbSC / HbAS 1.221 0.564 2.02
## HbSC / HbCC 0.776 0.238 1.74
##
## Point estimate displayed: median
## Results are back-transformed from the log odds ratio scale
## HPD interval probability: 0.95# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)
contrast(emm, method = list("HbAC / HbAA" = AC - AA,
"HbAS / HbAA" = AS - AA,
"HbCC / HbAA" = CC - AA,
"HbSC / HbAA" = SC - AA),
adjust = "bonferroni",
type="response",
infer = c(TRUE, TRUE)
)## contrast odds.ratio lower.HPD upper.HPD
## HbAC / HbAA 0.903 0.586 1.30
## HbAS / HbAA 0.783 0.499 1.10
## HbCC / HbAA 1.230 0.417 2.27
## HbSC / HbAA 0.955 0.466 1.60
##
## Point estimate displayed: median
## Results are back-transformed from the log odds ratio scale
## HPD interval probability: 0.95confint(emm, type="response", adjust="bonferroni")2.1.2 IgG+ asexual iRBC
The binomial model was not applicable.
model <- glm(
cbind(raw_igGpos_asexual,
raw_igGneg_asexual) ~ HBB.genotype,
family = binomial, data = df)
summary(model)##
## Call:
## glm(formula = cbind(raw_igGpos_asexual, raw_igGneg_asexual) ~
## HBB.genotype, family = binomial, data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.29328 0.01178 -109.789 < 2e-16 ***
## HBB.genotypeHbAC 0.01222 0.01680 0.727 0.467
## HBB.genotypeHbAS -0.27538 0.01688 -16.310 < 2e-16 ***
## HBB.genotypeHbCC 0.23987 0.03097 7.745 9.54e-15 ***
## HBB.genotypeHbSC -0.10229 0.02549 -4.013 6.00e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 7977.4 on 50 degrees of freedom
## Residual deviance: 7482.8 on 46 degrees of freedom
## AIC: 7892.4
##
## Number of Fisher Scoring iterations: 4# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)
# Test for overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 7.2991, p-value < 2.2e-16
## alternative hypothesis: two.sided# Test for zero-inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)
The beta-binomial GLM was used as the final model. The data was modelled using generalized linear model using the beta-binomial family and a single dispersion parameter to account for overdispersion.
model <- glmmTMB(
cbind(raw_igGpos_asexual,
raw_igGneg_asexual) ~ HBB.genotype,
family = betabinomial(link = "logit"),
data = df
)
summary(model)## Family: betabinomial ( logit )
## Formula: cbind(raw_igGpos_asexual, raw_igGneg_asexual) ~ HBB.genotype
## Data: df
##
## AIC BIC logLik -2*log(L) df.resid
## 725.7 737.3 -356.8 713.7 45
##
##
## Dispersion parameter for betabinomial family (): 16.4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.27711 0.15202 -8.401 <2e-16 ***
## HBB.genotypeHbAC -0.09105 0.21565 -0.422 0.673
## HBB.genotypeHbAS -0.30215 0.21096 -1.432 0.152
## HBB.genotypeHbCC 0.29470 0.40383 0.730 0.466
## HBB.genotypeHbSC -0.09908 0.32566 -0.304 0.761
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)
# Test for overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.83947, p-value = 0.48
## alternative hypothesis: two.sided# Test for zero-inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)
The estimated odds ratios (using Wald Z-tests) for the different genotype contrasts are as follows:
# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(model, ~ HBB.genotype)
# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## HbAC / HbAA 0.913 0.197 Inf 0.507 1.64 1 -0.422 0.9934
## HbAS / HbAA 0.739 0.156 Inf 0.416 1.31 1 -1.432 0.6067
## HbAS / HbAC 0.810 0.173 Inf 0.452 1.45 1 -0.989 0.8603
## HbCC / HbAA 1.343 0.542 Inf 0.446 4.04 1 0.730 0.9496
## HbCC / HbAC 1.471 0.596 Inf 0.487 4.44 1 0.952 0.8763
## HbCC / HbAS 1.816 0.732 Inf 0.605 5.45 1 1.482 0.5743
## HbSC / HbAA 0.906 0.295 Inf 0.373 2.20 1 -0.304 0.9981
## HbSC / HbAC 0.992 0.325 Inf 0.406 2.42 1 -0.025 1.0000
## HbSC / HbAS 1.225 0.397 Inf 0.506 2.97 1 0.627 0.9709
## HbSC / HbCC 0.675 0.319 Inf 0.186 2.45 1 -0.833 0.9206
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 5 estimates
## Intervals are back-transformed from the log odds ratio scale
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)
contrast(emm, method = list("HbAC / HbAA" = AC - AA,
"HbAS / HbAA" = AS - AA,
"HbCC / HbAA" = CC - AA,
"HbSC / HbAA" = SC - AA),
adjust = "bonferroni",
type="response",
infer = c(TRUE, TRUE)
)## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## HbAC / HbAA 0.913 0.197 Inf 0.533 1.56 1 -0.422 1.0000
## HbAS / HbAA 0.739 0.156 Inf 0.436 1.25 1 -1.432 0.6083
## HbCC / HbAA 1.343 0.542 Inf 0.490 3.68 1 0.730 1.0000
## HbSC / HbAA 0.906 0.295 Inf 0.402 2.04 1 -0.304 1.0000
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 4 estimates
## Intervals are back-transformed from the log odds ratio scale
## P value adjustment: bonferroni method for 4 tests
## Tests are performed on the log odds ratio scaleThe above table not only shows the estimated odds ratio and its
corresponding p-value, but also the confidence interval of this odds
ratio (= asymp.LCL and asymp.UCL).
These p-values are adjusted for multiple testing using the Bonferroni method.
The estimates are already transformed from the log-scale, so they can be directly interpreted as odds ratios (OR).
The estimated SCR per group is, adjusted for multiple testing using the Bonferroni method:
confint(emm, type="response", adjust="bonferroni")2.1.2.1 Bayesian estimates
brm_model <- brm(
bf(
raw_igGpos_asexual | trials(raw_igGpos_asexual + raw_igGneg_asexual) ~
HBB.genotype
),
family = beta_binomial(link = "logit"),
data = df,
iter = 5000,
chains = 4,
# control = list(adapt_delta = 0.95),
seed = 42,
cores = 8
)
## Compiling Stan program...
## Start samplingsummary(brm_model)## Family: beta_binomial
## Links: mu = logit
## Formula: raw_igGpos_asexual | trials(raw_igGpos_asexual + raw_igGneg_asexual) ~ HBB.genotype
## Data: df (Number of observations: 51)
## Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
## total post-warmup draws = 10000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.27 0.16 -1.60 -0.96 1.00 8837 6807
## HBB.genotypeHbAC -0.09 0.23 -0.53 0.36 1.00 8943 7844
## HBB.genotypeHbAS -0.30 0.22 -0.74 0.13 1.00 8743 7089
## HBB.genotypeHbCC 0.23 0.48 -0.76 1.09 1.00 9564 6469
## HBB.genotypeHbSC -0.14 0.35 -0.87 0.54 1.00 9960 7026
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi 15.07 3.05 9.63 21.51 1.00 10702 7405
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(brm_model, ~ HBB.genotype)
# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))## contrast odds.ratio lower.HPD upper.HPD
## HbAC / HbAA 0.915 0.545 1.37
## HbAS / HbAA 0.741 0.455 1.10
## HbAS / HbAC 0.807 0.473 1.19
## HbCC / HbAA 1.296 0.340 2.71
## HbCC / HbAC 1.416 0.385 2.96
## HbCC / HbAS 1.749 0.451 3.61
## HbSC / HbAA 0.888 0.370 1.59
## HbSC / HbAC 0.967 0.413 1.76
## HbSC / HbAS 1.192 0.511 2.14
## HbSC / HbCC 0.683 0.151 1.76
##
## Point estimate displayed: median
## Results are back-transformed from the log odds ratio scale
## HPD interval probability: 0.95# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)
contrast(emm, method = list("HbAC / HbAA" = AC - AA,
"HbAS / HbAA" = AS - AA,
"HbCC / HbAA" = CC - AA,
"HbSC / HbAA" = SC - AA),
adjust = "bonferroni",
type="response",
infer = c(TRUE, TRUE)
)## contrast odds.ratio lower.HPD upper.HPD
## HbAC / HbAA 0.915 0.545 1.37
## HbAS / HbAA 0.741 0.455 1.10
## HbCC / HbAA 1.296 0.340 2.71
## HbSC / HbAA 0.888 0.370 1.59
##
## Point estimate displayed: median
## Results are back-transformed from the log odds ratio scale
## HPD interval probability: 0.95confint(emm, type="response", adjust="bonferroni")2.1.3 IgG+ sexual iRBC
The binomial model was not applicable.
model <- glm(
cbind(raw_igGpos_sexual,
raw_igGneg_sexual) ~ HBB.genotype,
family = binomial, data = df)
summary(model)##
## Call:
## glm(formula = cbind(raw_igGpos_sexual, raw_igGneg_sexual) ~ HBB.genotype,
## family = binomial, data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.72119 0.03966 -93.816 < 2e-16 ***
## HBB.genotypeHbAC -0.03870 0.05682 -0.681 0.49589
## HBB.genotypeHbAS -0.15682 0.05613 -2.794 0.00521 **
## HBB.genotypeHbCC 0.22857 0.10188 2.243 0.02487 *
## HBB.genotypeHbSC 0.02220 0.08398 0.264 0.79154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 444.12 on 50 degrees of freedom
## Residual deviance: 426.00 on 46 degrees of freedom
## AIC: 715.45
##
## Number of Fisher Scoring iterations: 4# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)
# Test for overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 4.6992, p-value < 2.2e-16
## alternative hypothesis: two.sided# Test for zero-inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)
The beta-binomial GLM was used as the final model. The data was modelled using generalized linear model using the beta-binomial family and a single dispersion parameter to account for overdispersion.
model <- glmmTMB(
cbind(raw_igGpos_sexual,
raw_igGneg_sexual) ~ HBB.genotype,
family = betabinomial(link = "logit"),
data = df
)
summary(model)## Family: betabinomial ( logit )
## Formula: cbind(raw_igGpos_sexual, raw_igGneg_sexual) ~ HBB.genotype
## Data: df
##
## AIC BIC logLik -2*log(L) df.resid
## 454.8 466.4 -221.4 442.8 45
##
##
## Dispersion parameter for betabinomial family (): 235
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.69892 0.11681 -31.67 <2e-16 ***
## HBB.genotypeHbAC -0.11348 0.16838 -0.67 0.500
## HBB.genotypeHbAS -0.15115 0.16111 -0.94 0.348
## HBB.genotypeHbCC 0.27159 0.29402 0.92 0.356
## HBB.genotypeHbSC 0.03024 0.24244 0.12 0.901
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)
# Test for overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.78507, p-value = 0.416
## alternative hypothesis: two.sided# Test for zero-inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)
The estimated odds ratios (using Wald Z-tests) for the different genotype contrasts are as follows:
# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(model, ~ HBB.genotype)
# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## HbAC / HbAA 0.893 0.150 Inf 0.564 1.41 1 -0.674 0.9621
## HbAS / HbAA 0.860 0.139 Inf 0.554 1.33 1 -0.938 0.8820
## HbAS / HbAC 0.963 0.160 Inf 0.612 1.52 1 -0.227 0.9994
## HbCC / HbAA 1.312 0.386 Inf 0.588 2.93 1 0.924 0.8878
## HbCC / HbAC 1.470 0.436 Inf 0.654 3.30 1 1.297 0.6929
## HbCC / HbAS 1.526 0.447 Inf 0.687 3.39 1 1.444 0.5992
## HbSC / HbAA 1.031 0.250 Inf 0.532 2.00 1 0.125 0.9999
## HbSC / HbAC 1.155 0.284 Inf 0.590 2.26 1 0.585 0.9774
## HbSC / HbAS 1.199 0.289 Inf 0.621 2.31 1 0.753 0.9438
## HbSC / HbCC 0.786 0.270 Inf 0.307 2.01 1 -0.701 0.9563
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 5 estimates
## Intervals are back-transformed from the log odds ratio scale
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)
contrast(emm, method = list("HbAC / HbAA" = AC - AA,
"HbAS / HbAA" = AS - AA,
"HbCC / HbAA" = CC - AA,
"HbSC / HbAA" = SC - AA),
adjust = "bonferroni",
type="response",
infer = c(TRUE, TRUE)
)## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## HbAC / HbAA 0.893 0.150 Inf 0.586 1.36 1 -0.674 1.0000
## HbAS / HbAA 0.860 0.139 Inf 0.575 1.29 1 -0.938 1.0000
## HbCC / HbAA 1.312 0.386 Inf 0.630 2.73 1 0.924 1.0000
## HbSC / HbAA 1.031 0.250 Inf 0.563 1.89 1 0.125 1.0000
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 4 estimates
## Intervals are back-transformed from the log odds ratio scale
## P value adjustment: bonferroni method for 4 tests
## Tests are performed on the log odds ratio scaleThe above table not only shows the estimated odds ratio and its
corresponding p-value, but also the confidence interval of this odds
ratio (= asymp.LCL and asymp.UCL).
These p-values are adjusted for multiple testing using the Bonferroni method.
The estimates are already transformed from the log-scale, so they can be directly interpreted as odds ratios (OR).
The estimated SCR per group is, adjusted for multiple testing using the Bonferroni method:
confint(emm, type="response", adjust="bonferroni")2.1.3.1 Bayesian estimates
brm_model <- brm(
bf(
raw_igGpos_sexual | trials(raw_igGpos_sexual + raw_igGneg_sexual) ~
HBB.genotype
),
family = beta_binomial(link = "logit"),
data = df,
iter = 5000,
chains = 4,
# control = list(adapt_delta = 0.95),
seed = 42,
cores = 8
)
## Compiling Stan program...
## Start samplingsummary(brm_model)## Family: beta_binomial
## Links: mu = logit
## Formula: raw_igGpos_sexual | trials(raw_igGpos_sexual + raw_igGneg_sexual) ~ HBB.genotype
## Data: df (Number of observations: 51)
## Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
## total post-warmup draws = 10000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -3.69 0.13 -3.95 -3.44 1.00 8685 7215
## HBB.genotypeHbAC -0.12 0.19 -0.49 0.25 1.00 9199 7839
## HBB.genotypeHbAS -0.15 0.18 -0.49 0.20 1.00 9227 7465
## HBB.genotypeHbCC 0.20 0.35 -0.59 0.81 1.00 9808 6006
## HBB.genotypeHbSC -0.00 0.27 -0.59 0.50 1.00 10082 7296
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi 192.71 45.35 115.79 291.89 1.00 11297 7268
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(brm_model, ~ HBB.genotype)
# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))## contrast odds.ratio lower.HPD upper.HPD
## HbAC / HbAA 0.889 0.585 1.24
## HbAS / HbAA 0.862 0.597 1.20
## HbAS / HbAC 0.968 0.659 1.38
## HbCC / HbAA 1.249 0.486 2.14
## HbCC / HbAC 1.403 0.516 2.41
## HbCC / HbAS 1.450 0.580 2.49
## HbSC / HbAA 1.010 0.517 1.59
## HbSC / HbAC 1.137 0.585 1.83
## HbSC / HbAS 1.173 0.617 1.84
## HbSC / HbCC 0.811 0.292 1.68
##
## Point estimate displayed: median
## Results are back-transformed from the log odds ratio scale
## HPD interval probability: 0.95# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)
contrast(emm, method = list("HbAC / HbAA" = AC - AA,
"HbAS / HbAA" = AS - AA,
"HbCC / HbAA" = CC - AA,
"HbSC / HbAA" = SC - AA),
adjust = "bonferroni",
type="response",
infer = c(TRUE, TRUE)
)## contrast odds.ratio lower.HPD upper.HPD
## HbAC / HbAA 0.889 0.585 1.24
## HbAS / HbAA 0.862 0.597 1.20
## HbCC / HbAA 1.249 0.486 2.14
## HbSC / HbAA 1.010 0.517 1.59
##
## Point estimate displayed: median
## Results are back-transformed from the log odds ratio scale
## HPD interval probability: 0.95confint(emm, type="response", adjust="bonferroni")2.2 IgM
2.2.1 IgM+ iRBC
Proportion of IgM positive versus IgM negative infected red blood cells.
A binomial GLM is not applicable because of violated model assumptions.
model <- glm(
cbind(raw_igMpos_iRBC,
raw_igMneg_iRBC) ~ HBB.genotype,
family = binomial, data = df)
summary(model)##
## Call:
## glm(formula = cbind(raw_igMpos_iRBC, raw_igMneg_iRBC) ~ HBB.genotype,
## family = binomial, data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.355457 0.026214 -128.002 < 2e-16 ***
## HBB.genotypeHbAC 0.078928 0.036758 2.147 0.031773 *
## HBB.genotypeHbAS 0.004638 0.035976 0.129 0.897430
## HBB.genotypeHbCC -0.270520 0.081797 -3.307 0.000942 ***
## HBB.genotypeHbSC -0.610926 0.070615 -8.651 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2713.6 on 50 degrees of freedom
## Residual deviance: 2587.5 on 46 degrees of freedom
## AIC: 2911.4
##
## Number of Fisher Scoring iterations: 5# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)
# Test for overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 19.074, p-value < 2.2e-16
## alternative hypothesis: two.sided# Test for zero-inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)
The beta-binomial GLM was used as the final model. The data was modelled using generalized linear model using the beta-binomial family and a single dispersion parameter to account for overdispersion.
model <- glmmTMB(
cbind(raw_igMpos_iRBC,
raw_igMneg_iRBC) ~ HBB.genotype,
family = betabinomial(link = "logit"),
data = df
)
summary(model)## Family: betabinomial ( logit )
## Formula: cbind(raw_igMpos_iRBC, raw_igMneg_iRBC) ~ HBB.genotype
## Data: df
##
## AIC BIC logLik -2*log(L) df.resid
## 569.0 580.6 -278.5 557.0 45
##
##
## Dispersion parameter for betabinomial family (): 66.2
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.32915 0.16831 -19.780 <2e-16 ***
## HBB.genotypeHbAC -0.06698 0.23416 -0.286 0.775
## HBB.genotypeHbAS -0.03856 0.22235 -0.173 0.862
## HBB.genotypeHbCC -0.06662 0.47091 -0.141 0.888
## HBB.genotypeHbSC -0.43597 0.39062 -1.116 0.264
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)
# Test for overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.4807, p-value = 0.14
## alternative hypothesis: two.sided# Test for zero-inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0, p-value = 1
## alternative hypothesis: two.sided# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)
The estimated odds ratios (using Wald Z-tests) for the different genotype contrasts are as follows:
# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(model, ~ HBB.genotype)
# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## HbAC / HbAA 0.935 0.219 Inf 0.494 1.77 1 -0.286 0.9985
## HbAS / HbAA 0.962 0.214 Inf 0.525 1.76 1 -0.173 0.9998
## HbAS / HbAC 1.029 0.232 Inf 0.556 1.90 1 0.126 0.9999
## HbCC / HbAA 0.936 0.441 Inf 0.259 3.38 1 -0.141 0.9999
## HbCC / HbAC 1.000 0.472 Inf 0.276 3.63 1 0.001 1.0000
## HbCC / HbAS 0.972 0.454 Inf 0.272 3.47 1 -0.060 1.0000
## HbSC / HbAA 0.647 0.253 Inf 0.223 1.88 1 -1.116 0.7981
## HbSC / HbAC 0.691 0.271 Inf 0.237 2.02 1 -0.941 0.8810
## HbSC / HbAS 0.672 0.259 Inf 0.235 1.92 1 -1.031 0.8409
## HbSC / HbCC 0.691 0.391 Inf 0.147 3.24 1 -0.652 0.9663
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 5 estimates
## Intervals are back-transformed from the log odds ratio scale
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)
contrast(emm, method = list("HbAC / HbAA" = AC - AA,
"HbAS / HbAA" = AS - AA,
"HbCC / HbAA" = CC - AA,
"HbSC / HbAA" = SC - AA),
adjust = "bonferroni",
type="response",
infer = c(TRUE, TRUE)
)## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## HbAC / HbAA 0.935 0.219 Inf 0.521 1.68 1 -0.286 1.0000
## HbAS / HbAA 0.962 0.214 Inf 0.552 1.68 1 -0.173 1.0000
## HbCC / HbAA 0.936 0.441 Inf 0.289 3.03 1 -0.141 1.0000
## HbSC / HbAA 0.647 0.253 Inf 0.244 1.72 1 -1.116 1.0000
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 4 estimates
## Intervals are back-transformed from the log odds ratio scale
## P value adjustment: bonferroni method for 4 tests
## Tests are performed on the log odds ratio scaleThe above table not only shows the estimated odds ratio and its
corresponding p-value, but also the confidence interval of this odds
ratio (= asymp.LCL and asymp.UCL).
These p-values are adjusted for multiple testing using the Bonferroni method.
The estimates are already transformed from the log-scale, so they can be directly interpreted as odds ratios (OR).
The estimated SCR per group is, adjusted for multiple testing using the Bonferroni method:
confint(emm, type="response", adjust="bonferroni")2.2.1.1 Bayesian estimates
brm_model <- brm(
bf(
raw_igMpos_iRBC | trials(raw_igMpos_iRBC + raw_igMneg_iRBC) ~
HBB.genotype
),
family = beta_binomial(link = "logit"),
data = df,
iter = 5000,
chains = 4,
# control = list(adapt_delta = 0.95),
seed = 42,
cores = 8
)
## Compiling Stan program...
## Start samplingsummary(brm_model)## Family: beta_binomial
## Links: mu = logit
## Formula: raw_igMpos_iRBC | trials(raw_igMpos_iRBC + raw_igMneg_iRBC) ~ HBB.genotype
## Data: df (Number of observations: 51)
## Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
## total post-warmup draws = 10000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -3.32 0.18 -3.68 -2.98 1.00 8818 7200
## HBB.genotypeHbAC -0.06 0.25 -0.54 0.43 1.00 9089 7851
## HBB.genotypeHbAS -0.03 0.23 -0.48 0.44 1.00 9147 7077
## HBB.genotypeHbCC -0.23 0.57 -1.53 0.72 1.00 9201 5612
## HBB.genotypeHbSC -0.51 0.44 -1.47 0.28 1.00 9252 6355
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi 60.07 12.87 37.31 87.21 1.00 9920 6906
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(brm_model, ~ HBB.genotype)
# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))## contrast odds.ratio lower.HPD upper.HPD
## HbAC / HbAA 0.939 0.5227 1.44
## HbAS / HbAA 0.971 0.5785 1.48
## HbAS / HbAC 1.031 0.5993 1.55
## HbCC / HbAA 0.851 0.1254 1.85
## HbCC / HbAC 0.895 0.1284 1.97
## HbCC / HbAS 0.868 0.1642 1.92
## HbSC / HbAA 0.622 0.1656 1.21
## HbSC / HbAC 0.662 0.1874 1.29
## HbSC / HbAS 0.640 0.1845 1.23
## HbSC / HbCC 0.737 0.0844 2.45
##
## Point estimate displayed: median
## Results are back-transformed from the log odds ratio scale
## HPD interval probability: 0.95# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)
contrast(emm, method = list("HbAC / HbAA" = AC - AA,
"HbAS / HbAA" = AS - AA,
"HbCC / HbAA" = CC - AA,
"HbSC / HbAA" = SC - AA),
adjust = "bonferroni",
type="response",
infer = c(TRUE, TRUE)
)## contrast odds.ratio lower.HPD upper.HPD
## HbAC / HbAA 0.939 0.523 1.44
## HbAS / HbAA 0.971 0.578 1.48
## HbCC / HbAA 0.851 0.125 1.85
## HbSC / HbAA 0.622 0.166 1.21
##
## Point estimate displayed: median
## Results are back-transformed from the log odds ratio scale
## HPD interval probability: 0.95confint(emm, type="response", adjust="bonferroni")2.2.2 IgM+ asexual iRBC
The binomial model was not applicable.
model <- glm(
cbind(raw_igMpos_asexual,
raw_igMneg_asexual) ~ HBB.genotype,
family = binomial, data = df)
summary(model)##
## Call:
## glm(formula = cbind(raw_igMpos_asexual, raw_igMneg_asexual) ~
## HBB.genotype, family = binomial, data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.372290 0.026616 -126.700 < 2e-16 ***
## HBB.genotypeHbAC 0.068161 0.037417 1.822 0.06851 .
## HBB.genotypeHbAS -0.001644 0.036614 -0.045 0.96419
## HBB.genotypeHbCC -0.272719 0.083005 -3.286 0.00102 **
## HBB.genotypeHbSC -0.590718 0.071156 -8.302 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2613.4 on 50 degrees of freedom
## Residual deviance: 2500.1 on 46 degrees of freedom
## AIC: 2822.2
##
## Number of Fisher Scoring iterations: 5# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)
# Test for overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 19.528, p-value < 2.2e-16
## alternative hypothesis: two.sided# Test for zero-inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)
The beta-binomial GLM was used as the final model. The data was modelled using generalized linear model using the beta-binomial family and a single dispersion parameter to account for overdispersion.
model <- glmmTMB(
cbind(raw_igMpos_asexual,
raw_igMneg_asexual) ~ HBB.genotype,
family = betabinomial(link = "logit"),
data = df
)
summary(model)## Family: betabinomial ( logit )
## Formula: cbind(raw_igMpos_asexual, raw_igMneg_asexual) ~ HBB.genotype
## Data: df
##
## AIC BIC logLik -2*log(L) df.resid
## 565.0 576.6 -276.5 553.0 45
##
##
## Dispersion parameter for betabinomial family (): 68.2
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.34908 0.16756 -19.988 <2e-16 ***
## HBB.genotypeHbAC -0.07570 0.23357 -0.324 0.746
## HBB.genotypeHbAS -0.03852 0.22143 -0.174 0.862
## HBB.genotypeHbCC -0.06949 0.46938 -0.148 0.882
## HBB.genotypeHbSC -0.42438 0.38807 -1.094 0.274
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)
# Test for overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.533, p-value = 0.08
## alternative hypothesis: two.sided# Test for zero-inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0, p-value = 1
## alternative hypothesis: two.sided# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)
The estimated odds ratios (using Wald Z-tests) for the different genotype contrasts are as follows:
# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(model, ~ HBB.genotype)
# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## HbAC / HbAA 0.927 0.217 Inf 0.490 1.75 1 -0.324 0.9976
## HbAS / HbAA 0.962 0.213 Inf 0.526 1.76 1 -0.174 0.9998
## HbAS / HbAC 1.038 0.233 Inf 0.562 1.92 1 0.165 0.9998
## HbCC / HbAA 0.933 0.438 Inf 0.259 3.36 1 -0.148 0.9999
## HbCC / HbAC 1.006 0.474 Inf 0.278 3.64 1 0.013 1.0000
## HbCC / HbAS 0.969 0.451 Inf 0.273 3.45 1 -0.067 1.0000
## HbSC / HbAA 0.654 0.254 Inf 0.227 1.89 1 -1.094 0.8100
## HbSC / HbAC 0.706 0.275 Inf 0.244 2.04 1 -0.894 0.8991
## HbSC / HbAS 0.680 0.260 Inf 0.239 1.93 1 -1.008 0.8518
## HbSC / HbCC 0.701 0.395 Inf 0.151 3.26 1 -0.630 0.9703
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 5 estimates
## Intervals are back-transformed from the log odds ratio scale
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)
contrast(emm, method = list("HbAC / HbAA" = AC - AA,
"HbAS / HbAA" = AS - AA,
"HbCC / HbAA" = CC - AA,
"HbSC / HbAA" = SC - AA),
adjust = "bonferroni",
type="response",
infer = c(TRUE, TRUE)
)## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## HbAC / HbAA 0.927 0.217 Inf 0.517 1.66 1 -0.324 1.0000
## HbAS / HbAA 0.962 0.213 Inf 0.553 1.67 1 -0.174 1.0000
## HbCC / HbAA 0.933 0.438 Inf 0.289 3.01 1 -0.148 1.0000
## HbSC / HbAA 0.654 0.254 Inf 0.248 1.72 1 -1.094 1.0000
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 4 estimates
## Intervals are back-transformed from the log odds ratio scale
## P value adjustment: bonferroni method for 4 tests
## Tests are performed on the log odds ratio scaleThe above table not only shows the estimated odds ratio and its
corresponding p-value, but also the confidence interval of this odds
ratio (= asymp.LCL and asymp.UCL).
These p-values are adjusted for multiple testing using the Bonferroni method.
The estimates are already transformed from the log-scale, so they can be directly interpreted as odds ratios (OR).
The estimated SCR per group is, adjusted for multiple testing using the Bonferroni method:
confint(emm, type="response", adjust="bonferroni")2.2.2.1 Bayesian estimates
brm_model <- brm(
bf(
raw_igMpos_asexual | trials(raw_igMpos_asexual + raw_igMneg_asexual) ~
HBB.genotype
),
family = beta_binomial(link = "logit"),
data = df,
iter = 5000,
chains = 4,
# control = list(adapt_delta = 0.95),
seed = 42,
cores = 8
)
## Compiling Stan program...
## Start samplingsummary(brm_model)## Family: beta_binomial
## Links: mu = logit
## Formula: raw_igMpos_asexual | trials(raw_igMpos_asexual + raw_igMneg_asexual) ~ HBB.genotype
## Data: df (Number of observations: 51)
## Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
## total post-warmup draws = 10000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -3.34 0.18 -3.71 -3.00 1.00 8254 6656
## HBB.genotypeHbAC -0.07 0.25 -0.56 0.42 1.00 8688 7179
## HBB.genotypeHbAS -0.03 0.23 -0.47 0.44 1.00 8979 7311
## HBB.genotypeHbCC -0.24 0.58 -1.60 0.70 1.00 9310 5117
## HBB.genotypeHbSC -0.49 0.43 -1.44 0.28 1.00 9819 5964
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi 61.74 13.42 38.01 91.19 1.00 9051 6501
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(brm_model, ~ HBB.genotype)
# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))## contrast odds.ratio lower.HPD upper.HPD
## HbAC / HbAA 0.932 0.5367 1.45
## HbAS / HbAA 0.968 0.5748 1.46
## HbAS / HbAC 1.043 0.6174 1.57
## HbCC / HbAA 0.842 0.1035 1.79
## HbCC / HbAC 0.899 0.1303 1.97
## HbCC / HbAS 0.866 0.1012 1.84
## HbSC / HbAA 0.635 0.1707 1.20
## HbSC / HbAC 0.680 0.2016 1.30
## HbSC / HbAS 0.655 0.1688 1.20
## HbSC / HbCC 0.759 0.0971 2.45
##
## Point estimate displayed: median
## Results are back-transformed from the log odds ratio scale
## HPD interval probability: 0.95# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)
contrast(emm, method = list("HbAC / HbAA" = AC - AA,
"HbAS / HbAA" = AS - AA,
"HbCC / HbAA" = CC - AA,
"HbSC / HbAA" = SC - AA),
adjust = "bonferroni",
type="response",
infer = c(TRUE, TRUE)
)## contrast odds.ratio lower.HPD upper.HPD
## HbAC / HbAA 0.932 0.537 1.45
## HbAS / HbAA 0.968 0.575 1.46
## HbCC / HbAA 0.842 0.103 1.79
## HbSC / HbAA 0.635 0.171 1.20
##
## Point estimate displayed: median
## Results are back-transformed from the log odds ratio scale
## HPD interval probability: 0.95confint(emm, type="response", adjust="bonferroni")2.2.3 IgM+ sexual iRBC
The binomial model was not applicable.
model <- glm(
cbind(raw_igMpos_sexual,
raw_igMneg_sexual) ~ HBB.genotype,
family = binomial, data = df)
summary(model)##
## Call:
## glm(formula = cbind(raw_igMpos_sexual, raw_igMneg_sexual) ~ HBB.genotype,
## family = binomial, data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.1787 0.2183 -32.884 <2e-16 ***
## HBB.genotypeHbAC 0.4958 0.2777 1.785 0.0742 .
## HBB.genotypeHbAS 0.1753 0.2867 0.611 0.5409
## HBB.genotypeHbCC -1.1068 1.0237 -1.081 0.2796
## HBB.genotypeHbSC -17.0827 1278.5349 -0.013 0.9893
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 166.65 on 50 degrees of freedom
## Residual deviance: 146.32 on 46 degrees of freedom
## AIC: 227.54
##
## Number of Fisher Scoring iterations: 16# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)
# Test for overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 2.9299, p-value < 2.2e-16
## alternative hypothesis: two.sided# Test for zero-inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 2.0402, p-value < 2.2e-16
## alternative hypothesis: two.sided# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)
The beta-binomial GLM was used as the final model. The data was modelled using generalized linear model using the beta-binomial family and a single dispersion parameter to account for overdispersion.
model <- glmmTMB(
cbind(raw_igMpos_sexual,
raw_igMneg_sexual) ~ HBB.genotype,
family = betabinomial(link = "logit"),
data = df
)
summary(model)## Family: betabinomial ( logit )
## Formula: cbind(raw_igMpos_sexual, raw_igMneg_sexual) ~ HBB.genotype
## Data: df
##
## AIC BIC logLik -2*log(L) df.resid
## 180.0 191.5 -84.0 168.0 45
##
##
## Dispersion parameter for betabinomial family (): 573
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.007e+00 3.688e-01 -18.999 <2e-16 ***
## HBB.genotypeHbAC 1.406e-01 4.868e-01 0.289 0.773
## HBB.genotypeHbAS 2.055e-02 4.631e-01 0.044 0.965
## HBB.genotypeHbCC -4.401e-01 1.057e+00 -0.416 0.677
## HBB.genotypeHbSC -1.974e+01 1.107e+04 -0.002 0.999
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)
# Test for overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.80213, p-value = 0.83
## alternative hypothesis: two.sided# Test for zero-inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0213, p-value = 0.92
## alternative hypothesis: two.sided# plot residuals per group
plotResiduals(sim_res, df$HBB.genotype)
The estimated odds ratios (using Wald Z-tests) for the different genotype contrasts are as follows:
# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(model, ~ HBB.genotype)
# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## HbAC / HbAA 1.151 5.60e-01 Inf 0.3051 4 1 0.289 0.9985
## HbAS / HbAA 1.021 4.73e-01 Inf 0.2886 4 1 0.044 1.0000
## HbAS / HbAC 0.887 4.17e-01 Inf 0.2462 3 1 -0.256 0.9991
## HbCC / HbAA 0.644 6.81e-01 Inf 0.0360 12 1 -0.416 0.9937
## HbCC / HbAC 0.560 5.95e-01 Inf 0.0307 10 1 -0.546 0.9825
## HbCC / HbAS 0.631 6.63e-01 Inf 0.0358 11 1 -0.438 0.9924
## HbSC / HbAA 0.000 2.95e-05 Inf 0.0000 Inf 1 -0.002 1.0000
## HbSC / HbAC 0.000 2.56e-05 Inf 0.0000 Inf 1 -0.002 1.0000
## HbSC / HbAS 0.000 2.89e-05 Inf 0.0000 Inf 1 -0.002 1.0000
## HbSC / HbCC 0.000 4.57e-05 Inf 0.0000 Inf 1 -0.002 1.0000
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 5 estimates
## Intervals are back-transformed from the log odds ratio scale
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)
contrast(emm, method = list("HbAC / HbAA" = AC - AA,
"HbAS / HbAA" = AS - AA,
"HbCC / HbAA" = CC - AA,
"HbSC / HbAA" = SC - AA),
adjust = "bonferroni",
type="response",
infer = c(TRUE, TRUE)
)## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## HbAC / HbAA 1.151 5.60e-01 Inf 0.3412 4 1 0.289 1.0000
## HbAS / HbAA 1.021 4.73e-01 Inf 0.3210 3 1 0.044 1.0000
## HbCC / HbAA 0.644 6.81e-01 Inf 0.0459 9 1 -0.416 1.0000
## HbSC / HbAA 0.000 2.95e-05 Inf 0.0000 Inf 1 -0.002 1.0000
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 4 estimates
## Intervals are back-transformed from the log odds ratio scale
## P value adjustment: bonferroni method for 4 tests
## Tests are performed on the log odds ratio scaleThe above table not only shows the estimated odds ratio and its
corresponding p-value, but also the confidence interval of this odds
ratio (= asymp.LCL and asymp.UCL).
These p-values are adjusted for multiple testing using the Bonferroni method.
The estimates are already transformed from the log-scale, so they can be directly interpreted as odds ratios (OR).
The estimated SCR per group is, adjusted for multiple testing using the Bonferroni method:
confint(emm, type="response", adjust="bonferroni")2.2.3.1 Bayesian estimates
brm_model <- brm(
bf(
raw_igMpos_sexual | trials(raw_igMpos_sexual + raw_igMneg_sexual) ~
HBB.genotype
),
family = beta_binomial(link = "logit"),
data = df,
iter = 5000,
chains = 4,
# control = list(adapt_delta = 0.95),
seed = 42,
cores = 8
)
## Compiling Stan program...
## Start sampling
## Warning: There were 37 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problemssummary(brm_model)## Warning: There were 37 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Family: beta_binomial
## Links: mu = logit
## Formula: raw_igMpos_sexual | trials(raw_igMpos_sexual + raw_igMneg_sexual) ~ HBB.genotype
## Data: df (Number of observations: 51)
## Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
## total post-warmup draws = 10000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -6.72 0.41 -7.55 -5.94 1.00 3695 4338
## HBB.genotypeHbAC 0.04 0.51 -0.98 1.04 1.00 3673 3306
## HBB.genotypeHbAS -0.01 0.49 -0.97 0.98 1.00 3907 3915
## HBB.genotypeHbCC -0.85 1.29 -4.01 1.16 1.00 3809 3003
## HBB.genotypeHbSC -57.62 72.78 -270.09 -2.42 1.00 1057 711
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi 294.73 114.34 112.63 561.06 1.00 4388 4357
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Compute estimated marginal means for the HBB.genotype factor
emm <- emmeans(brm_model, ~ HBB.genotype)
# Perform all pairwise comparisons
pairs(emm, adjust="tukey", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))## contrast odds.ratio lower.HPD upper.HPD
## HbAC / HbAA 1.043 0.277451 2.4730
## HbAS / HbAA 0.984 0.262777 2.2540
## HbAS / HbAC 0.950 0.239007 2.2416
## HbCC / HbAA 0.520 0.000906 2.4482
## HbCC / HbAC 0.495 0.000608 2.3643
## HbCC / HbAS 0.520 0.000762 2.4449
## HbSC / HbAA 0.000 0.000000 0.0272
## HbSC / HbAC 0.000 0.000000 0.0264
## HbSC / HbAS 0.000 0.000000 0.0277
## HbSC / HbCC 0.000 0.000000 0.0702
##
## Point estimate displayed: median
## Results are back-transformed from the log odds ratio scale
## HPD interval probability: 0.95# perform only comparisons against AA
AA = c(1,0,0,0,0)
AC = c(0,1,0,0,0)
AS = c(0,0,1,0,0)
CC = c(0,0,0,1,0)
SC = c(0,0,0,0,1)
contrast(emm, method = list("HbAC / HbAA" = AC - AA,
"HbAS / HbAA" = AS - AA,
"HbCC / HbAA" = CC - AA,
"HbSC / HbAA" = SC - AA),
adjust = "bonferroni",
type="response",
infer = c(TRUE, TRUE)
)## contrast odds.ratio lower.HPD upper.HPD
## HbAC / HbAA 1.043 0.277451 2.4730
## HbAS / HbAA 0.984 0.262777 2.2540
## HbCC / HbAA 0.520 0.000906 2.4482
## HbSC / HbAA 0.000 0.000000 0.0272
##
## Point estimate displayed: median
## Results are back-transformed from the log odds ratio scale
## HPD interval probability: 0.95confint(emm, type="response", adjust="bonferroni")2.3 Pf infection
2.3.1 IgG+
Proportion of IgG positive versus IgG negative infected red blood cells for active vs non-active infections:
A binomial GLM is not applicable because of violated model assumptions.
model <- glm(
cbind(raw_igGpos_iRBC,
raw_igGneg_iRBC) ~ Pf_infection,
family = binomial, data = df)
summary(model)##
## Call:
## glm(formula = cbind(raw_igGpos_iRBC, raw_igGneg_iRBC) ~ Pf_infection,
## family = binomial, data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.96465 0.01127 -174.35 <2e-16 ***
## Pf_infectionYES 0.48053 0.01284 37.42 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 8242.3 on 50 degrees of freedom
## Residual deviance: 6752.8 on 49 degrees of freedom
## AIC: 7176.8
##
## Number of Fisher Scoring iterations: 4# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)
# Test for overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 4.6219, p-value < 2.2e-16
## alternative hypothesis: two.sided# Test for zero-inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided# plot residuals per group
plotResiduals(sim_res, df$Pf_infection)
The beta-binomial GLM was used as the final model. The data was modelled using generalized linear model using the beta-binomial family and a single dispersion parameter to account for overdispersion.
model <- glmmTMB(
cbind(raw_igGpos_iRBC,
raw_igGneg_iRBC) ~ Pf_infection,
family = betabinomial(link = "logit"),
data = df
)
summary(model)## Family: betabinomial ( logit )
## Formula: cbind(raw_igGpos_iRBC, raw_igGneg_iRBC) ~ Pf_infection
## Data: df
##
## AIC BIC logLik -2*log(L) df.resid
## 728.2 734.0 -361.1 722.2 48
##
##
## Dispersion parameter for betabinomial family (): 34.8
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.0615 0.1327 -15.532 < 2e-16 ***
## Pf_infectionYES 0.5944 0.1490 3.989 6.63e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)
# Test for overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.87804, p-value = 0.568
## alternative hypothesis: two.sided# Test for zero-inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided# plot residuals per group
plotResiduals(sim_res, df$Pf_infection)
The estimated odds ratios (using Wald Z-tests) for the different genotype contrasts are as follows:
# Compute estimated marginal means for the Pf_infection factor
emm <- emmeans(model, ~ Pf_infection)
# Perform all pairwise comparisons
pairs(emm, adjust="bonferroni", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## YES / NO 1.81 0.27 Inf 1.35 2.43 1 3.989 0.0001
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log odds ratio scale
## Tests are performed on the log odds ratio scaleThe above table not only shows the estimated odds ratio and its
corresponding p-value, but also the confidence interval of this odds
ratio (= asymp.LCL and asymp.UCL).
These p-values are adjusted for multiple testing using the Bonferroni method.
The estimates are already transformed from the log-scale, so they can be directly interpreted as odds ratios (OR).
The estimated SCR per group is, adjusted for multiple testing using the Bonferroni method:
confint(emm, type="response", adjust="bonferroni")2.3.1.1 Bayesian estimates
brm_model <- brm(
bf(
raw_igGpos_iRBC | trials(raw_igGpos_iRBC + raw_igGneg_iRBC) ~
Pf_infection
),
family = beta_binomial(link = "logit"),
data = df,
iter = 5000,
chains = 4,
# control = list(adapt_delta = 0.95),
seed = 42,
cores = 8
)
## Compiling Stan program...
## Start samplingsummary(brm_model)## Family: beta_binomial
## Links: mu = logit
## Formula: raw_igGpos_iRBC | trials(raw_igGpos_iRBC + raw_igGneg_iRBC) ~ Pf_infection
## Data: df (Number of observations: 51)
## Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
## total post-warmup draws = 10000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -2.06 0.14 -2.34 -1.80 1.00 6802 6563
## Pf_infectionYES 0.60 0.16 0.30 0.91 1.00 7860 6694
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi 33.18 6.71 21.18 47.73 1.00 7335 6405
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Compute estimated marginal means for the Pf_infection factor
emm <- emmeans(brm_model, ~ Pf_infection)
# Perform all pairwise comparisons
pairs(emm, adjust="bonferroni", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))## contrast odds.ratio lower.HPD upper.HPD
## YES / NO 1.82 1.31 2.43
##
## Point estimate displayed: median
## Results are back-transformed from the log odds ratio scale
## HPD interval probability: 0.95confint(emm, type="response", adjust="bonferroni")2.3.1.2 Model comparison
model <- glmmTMB(
cbind(raw_igGpos_iRBC,
raw_igGneg_iRBC) ~ Pf_infection,
family = betabinomial(link = "logit"),
data = df
)
model_extended <- glmmTMB(
cbind(raw_igGpos_iRBC,
raw_igGneg_iRBC) ~ Pf_infection + Gender + Village + Age_group,
family = betabinomial(link = "logit"),
data = df
)
anova(model, model_extended)A more complex model is not preferred. Moreover, the directionality and magnitude of the effect remains consistent.
summary(model_extended)## Family: betabinomial ( logit )
## Formula:
## cbind(raw_igGpos_iRBC, raw_igGneg_iRBC) ~ Pf_infection + Gender +
## Village + Age_group
## Data: df
##
## AIC BIC logLik -2*log(L) df.resid
## 728.5 743.9 -356.2 712.5 43
##
##
## Dispersion parameter for betabinomial family (): 42.2
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.90038 0.15024 -12.649 < 2e-16 ***
## Pf_infectionYES 0.66309 0.14878 4.457 8.32e-06 ***
## Gender1 -0.03191 0.11609 -0.275 0.78339
## Village2 -0.33062 0.13737 -2.407 0.01609 *
## Village3 -0.23119 0.20133 -1.148 0.25083
## Village4 -0.52382 0.19835 -2.641 0.00827 **
## Age_group2 -0.03899 0.12991 -0.300 0.76409
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Simulate residuals
sim_res <- simulateResiduals(model_extended, n = 1000)
plot(sim_res)
# Test for overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84876, p-value = 0.462
## alternative hypothesis: two.sided# Test for zero-inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided# plot residuals per group
plotResiduals(sim_res, df$Pf_infection)
2.3.2 IgM+
Proportion of IgM positive versus IgM negative infected red blood cells for active vs non-active infections:
A binomial GLM is not applicable because of violated model assumptions.
model <- glm(
cbind(raw_igMpos_iRBC,
raw_igMneg_iRBC) ~ Pf_infection,
family = binomial, data = df)
summary(model)##
## Call:
## glm(formula = cbind(raw_igMpos_iRBC, raw_igMneg_iRBC) ~ Pf_infection,
## family = binomial, data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.87580 0.03309 -117.11 <2e-16 ***
## Pf_infectionYES 0.64817 0.03660 17.71 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2713.6 on 50 degrees of freedom
## Residual deviance: 2358.5 on 49 degrees of freedom
## AIC: 2676.5
##
## Number of Fisher Scoring iterations: 5# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)
# Test for overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 8.3282, p-value < 2.2e-16
## alternative hypothesis: two.sided# Test for zero-inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided# plot residuals per group
plotResiduals(sim_res, df$Pf_infection)
The beta-binomial GLM was used as the final model. The data was modelled using generalized linear model using the beta-binomial family and a single dispersion parameter to account for overdispersion.
model <- glmmTMB(
cbind(raw_igMpos_iRBC,
raw_igMneg_iRBC) ~ Pf_infection,
family = betabinomial(link = "logit"),
data = df
)
summary(model)## Family: betabinomial ( logit )
## Formula: cbind(raw_igMpos_iRBC, raw_igMneg_iRBC) ~ Pf_infection
## Data: df
##
## AIC BIC logLik -2*log(L) df.resid
## 560.1 565.9 -277.1 554.1 48
##
##
## Dispersion parameter for betabinomial family (): 70.5
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.6959 0.1830 -20.192 <2e-16 ***
## Pf_infectionYES 0.4089 0.2016 2.029 0.0425 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Simulate residuals
sim_res <- simulateResiduals(model, n = 1000)
plot(sim_res)
# Test for overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.5207, p-value = 0.126
## alternative hypothesis: two.sided# Test for zero-inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0, p-value = 1
## alternative hypothesis: two.sided# plot residuals per group
plotResiduals(sim_res, df$Pf_infection)
The estimated odds ratios (using Wald Z-tests) for the different genotype contrasts are as follows:
# Compute estimated marginal means for the Pf_infection factor
emm <- emmeans(model, ~ Pf_infection)
# Perform all pairwise comparisons
pairs(emm, adjust="bonferroni", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## YES / NO 1.51 0.303 Inf 1.01 2.23 1 2.029 0.0425
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log odds ratio scale
## Tests are performed on the log odds ratio scaleThe above table not only shows the estimated odds ratio and its
corresponding p-value, but also the confidence interval of this odds
ratio (= asymp.LCL and asymp.UCL).
These p-values are adjusted for multiple testing using the Bonferroni method.
The estimates are already transformed from the log-scale, so they can be directly interpreted as odds ratios (OR).
The estimated SCR per group is, adjusted for multiple testing using the Bonferroni method:
confint(emm, type="response", adjust="bonferroni")2.3.2.1 Bayesian estimates
brm_model <- brm(
bf(
raw_igMpos_iRBC | trials(raw_igMpos_iRBC + raw_igMneg_iRBC) ~
Pf_infection
),
family = beta_binomial(link = "logit"),
data = df,
iter = 5000,
chains = 4,
# control = list(adapt_delta = 0.95),
seed = 42,
cores = 8
)
## Compiling Stan program...
## Start samplingsummary(brm_model)## Family: beta_binomial
## Links: mu = logit
## Formula: raw_igMpos_iRBC | trials(raw_igMpos_iRBC + raw_igMneg_iRBC) ~ Pf_infection
## Data: df (Number of observations: 51)
## Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
## total post-warmup draws = 10000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -3.69 0.19 -4.07 -3.34 1.00 6274 6111
## Pf_infectionYES 0.42 0.21 0.03 0.83 1.00 7019 6525
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi 66.14 14.05 41.43 95.92 1.00 5906 5889
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Compute estimated marginal means for the Pf_infection factor
emm <- emmeans(brm_model, ~ Pf_infection)
# Perform all pairwise comparisons
pairs(emm, adjust="bonferroni", ratios=T, type="response", reverse=T, infer = c(TRUE, TRUE))## contrast odds.ratio lower.HPD upper.HPD
## YES / NO 1.51 0.963 2.2
##
## Point estimate displayed: median
## Results are back-transformed from the log odds ratio scale
## HPD interval probability: 0.95confint(emm, type="response", adjust="bonferroni")3 Software used
grateful::cite_packages(output = "table", out.dir = here(), cite.tidyverse = TRUE)